compare_datasets
Purpose
This markdown follows differential expression analyses from deseq2_analyses.Rmd. DEGs and signatures from pre-cystic P70 (‘39’) data and cystic P21 (‘56’) and P28 (‘19’) data are compared to each other and used for functional enirichment analysis
Data Sets - ’39 / Pre-Cystic P70: GSE149739
- ’19 / Cystic P21: GSE134719
- ’56 / Cystic P28: GSE69556
library(gprofiler2)
library(viridis)
library(stringr)
library(signatureSearch)
library(ggplot2)
library(dplyr)
library(here)
library(viridis)
library(rrvgo)
library(GOSemSim)
library(VennDiagram)
library(RColorBrewer)
library(ComplexHeatmap)
library(ExperimentHub)
source(here("code", "functions.R"))##Signature comparisons
Load in (top 100 human) data
Load in top 100 up and down genes used as signatures for signatureSearch
’39 precystic p70
hom_p70_UP <- read.csv(here("res", "deseq2_outputs", "p70_signature_UP.csv"), row.names = 1)
hom_p70_DOWN <- read.csv(here("res", "deseq2_outputs", "p70_signature_DOWN.csv"), row.names = 1)’19 cystic p28
hom_p28_UP <- read.csv(here("res", "deseq2_outputs", "p28_signature_UP.csv"), row.names = 1)
hom_p28_DOWN <- read.csv(here("res", "deseq2_outputs", "p28_signature_DOWN.csv"), row.names = 1)’56 cystic p21
hom_p21_UP <- read.csv(here("res", "deseq2_outputs", "p21_signature_UP.csv"), row.names = 1)
hom_p21_DOWN <- read.csv(here("res", "deseq2_outputs", "p21_signature_DOWN.csv"), row.names = 1)specific to pre-cystic
p70_unique_up <- dplyr::anti_join(hom_p70_UP, hom_p28_UP, by = "hgnc_symbol")
p70_unique_up <- dplyr::anti_join(p70_unique_up, hom_p21_UP, by = "hgnc_symbol")
p70_unique_down <- dplyr::anti_join(hom_p70_DOWN, hom_p28_DOWN, by = "hgnc_symbol")
p70_unique_down <- dplyr::anti_join(p70_unique_down, hom_p21_DOWN, by = "hgnc_symbol")up_p28_p21 <- intersect(hom_p28_UP$hgnc_symbol, hom_p21_UP$hgnc_symbol)
down_p28_p21 <- intersect(hom_p28_DOWN$hgnc_symbol, hom_p21_DOWN$hgnc_symbol)
up_all <- intersect(up_p28_p21, hom_p70_UP$hgnc_symbol)
down_all <- intersect(down_p28_p21, hom_p70_DOWN$hgnc_symbol)Venn diagram of signatures
setp70 <- c(hom_p70_UP$entrezgene_id, hom_p70_DOWN$entrezgene_id)
setp21 <- c(hom_p21_UP$entrezgene_id, hom_p21_DOWN$entrezgene_id)
setp28 <- c(hom_p28_UP$entrezgene_id, hom_p28_DOWN$entrezgene_id)
venn3(setp70, setp21, setp28, categories = c("Precystic P70" , "Cystic P21" , "Cystic P28"), filename = here("res/deseq2_outputs", "signatures_venndiagram.png"))## [1] 1FEA p70 precystic GSE1497p70 (top 100 human orthologs)
fea_p70 <- enrich(as.character(hom_p70_UP$entrezgene_id), as.character(hom_p70_DOWN$entrezgene_id), 11) Bubbleplot
bubbleplot(fea_p70$topcompared)
fea_p70$topcompared$term_name <- stringr::str_wrap(fea_p70$topcompared$term_name, width = 70)
bp_p70 <- bubbleplot(fea_p70$topcompared) + labs(title = "Top Enriched Terms for \nPrecystic P70 Signature") + scale_x_discrete(labels=c("downregulated" = "down", "upregulated" = "up")) + theme(text = element_text(size = 20, face = "bold"), plot.title = element_text(size = 20))
bp_p70
ggsave(bp_p70, filename = here("res", "fea", "bubbleplot_signature_p70.png"), height = 8
, width = 8.75)FEA p21 cystic GSE695p21 (top 100 human orthologs)
fea_p21 <- enrich(as.character(hom_p21_UP$entrezgene_id), as.character(hom_p21_DOWN$entrezgene_id), 9)fea_p21$topcompared$term_name <- stringr::str_trunc(fea_p21$topcompared$term_name, width = 55)
bp_p21 <- bubbleplot(fea_p21$topcompared) + labs(title = "Top Enriched Terms for \nCystic P21 Signature") + scale_x_discrete(labels=c("downregulated" = "down", "upregulated" = "up")) + theme(text = element_text(size = 20, face = "bold"), plot.title = element_text(size = 20))
bp_p21
ggsave(bp_p21, filename = here("res", "fea", "bubbleplot_signature_p21.png"), height = 9, width = 9.76)FEA p28 cystic GSE1347p28 (top 100 human orthologs)
fea_p28 <- enrich(as.character(hom_p28_UP$entrezgene_id), as.character(hom_p28_DOWN$entrezgene_id), 11)bubbleplot(fea_p28$topcompared)
fea_p28$topcompared$term_name <- stringr::str_trunc(fea_p28$topcompared$term_name, width = 42)
bp_p28 <- bubbleplot(fea_p28$topcompared) + labs(title = "Top Enriched Terms for \nCystic P28 Signature") + scale_x_discrete(labels=c("downregulated" = "down", "upregulated" = "up")) + theme(text = element_text(size = 20, face = "bold"), plot.title = element_text(size = 20))
bp_p28
ggsave(bp_p28, filename = here("res", "fea", "bubbleplot_signature_p28.png"), height = 8.7, width = 8.75)Save
##fea_p70
fea_p70_upterms <- apply(fea_p70$upregulated, 2, as.character)
fea_p70_downterms <- apply(fea_p70$downregulated, 2, as.character)
write.csv(rbind(fea_p70_upterms, fea_p70_downterms), file = here("res", "fea", "fea_p70sig.csv"))
##fea_p21
fea_p21_upterms <- apply(fea_p21$upregulated, 2, as.character)
fea_p21_downterms <- apply(fea_p21$downregulated, 2, as.character)
write.csv(rbind(fea_p21_upterms, fea_p21_downterms), file = here("res", "fea", "fea_p21sig.csv"))
##fea_p28
fea_p28_upterms <- apply(fea_p28$upregulated, 2, as.character)
fea_p28_downterms <- apply(fea_p28$downregulated, 2, as.character)
write.csv(rbind(fea_p28_upterms, fea_p28_downterms), file = here("res", "fea", "fea_p28sig.csv"))Term Overlaps
Term overlap
cystic_term_overlap_up <- merge(fea_p28$upregulated, fea_p21$upregulated, by = "term_name")
cystic_term_overlap_down <- merge(fea_p28$downregulated, fea_p21$downregulated, by = "term_name")
p70_p21_term_overlap_up <- merge(fea_p70$upregulated, fea_p21$upregulated, by = "term_name")
p70_p21_term_overlap_down <- merge(fea_p70$downregulated, fea_p21$downregulated, by = "term_name")
p70_p28_term_overlap_up <- merge(fea_p70$upregulated, fea_p28$upregulated, by = "term_name")
p70_p28_term_overlap_down <- merge(fea_p28$downregulated, fea_p28$downregulated, by = "term_name")
fea_overap_all_up <- merge(fea_p70$upregulated, cystic_term_overlap_up, by = "term_name")
fea_overap_all_down <- merge(fea_p70$downregulated, cystic_term_overlap_down, by = "term_name")Pathway Similarity
rrvgo, uses GO terms as input ##### P70
sim_p70 <- go_bp_sim(fea_p70, "org.Hs.eg.db", threshold = 0.80)## ## preparing gene to GO mapping data...## preparing IC data...## Warning in calculateSimMatrix(go_comb$term_id, orgdb = orgdb, ont = "BP", :
## Removed 1 terms that were not found in orgdb for BPsp <- scatterPlot2(sim_p70$simMatrix, sim_p70$reducedTerms, labelSize = 2.5) + labs(title = "Precystic P70 Signature Enrichment Overview") #+ geom_label_repel(max.overlaps = 20) #+ theme(text = element_text(size = 14), plot.title = element_text(size = 24))
sp
#tm <- treemapPlot(sim_p70$reducedTerms, title = "Enrichment for Precystic P70 Signature")
library(treemap)
tm2 <- treemapPlot2(sim_p70$reducedTerms, title = "Enrichment for Precystic P70 Signature") 
hm <- heatmapPlot(sim_p70$simMatrix, sim_p70$reducedTerms, annotateParent=TRUE, annotationLabel="parentTerm", fontsize=6)
#ggsave(here("res/fea", "precysticp70_sig_enrichment_scatter.pdf"), sp, height = 5, width = 7)
#ggsave(here("res/fea", "precysticp70_sig_enrichment_treemap.pdf"), tm, height = 4, width = 7)P21
sim_p21 <- go_bp_sim(fea_p21, "org.Hs.eg.db", threshold = 0.80)## preparing gene to GO mapping data...## preparing IC data...#options(ggrepel.max.overlaps = 50)
sp_p21 <- scatterPlot2(sim_p21$simMatrix, sim_p21$reducedTerms, labelSize = 2.5) + labs(title = "Cystic P21 Signature Enrichment Overview")
sp_p21
tm_p21 <- treemapPlot2(sim_p21$reducedTerms, title = "Enrichment for Cystic P21 Signature")
hm_p21 <- heatmapPlot(sim_p21$simMatrix, sim_p21$reducedTerms, annotateParent=TRUE, annotationLabel="parentTerm", fontsize=6)
ggsave(here("res/fea", "cysticp21_sig_enrichment_scatter.pdf"), sp_p21, height =5, width = 7)
#ggsave(here("res/fea", "precysticp70_sig_enrichment_treemap.pdf"), tm, height = 4, width = 7)P28
sim_p28 <- go_bp_sim(fea_p28, "org.Hs.eg.db", threshold = 0.80)## preparing gene to GO mapping data...## preparing IC data...#options(ggrepel.max.overlaps = 50)
sp_p28 <- scatterPlot2(sim_p28$simMatrix, sim_p28$reducedTerms, labelSize = 2.5) + labs(title = "Cystic P21 Signature Enrichment Overview")
sp_p28
tm_p28 <- treemapPlot2(sim_p28$reducedTerms, title = "Enrichment for Cystic P28 Signature")
hm_p28 <- heatmapPlot(sim_p28$simMatrix, sim_p28$reducedTerms, annotateParent=TRUE, annotationLabel="parentTerm", fontsize=6)
ggsave(here("res/fea", "cysticp28_sig_enrichment_scatter.pdf"), sp_p28, height =5, width = 7)DEG Comparisons
Using ALL mouse DEGs (0.05 alpha and -2 and 2 LFC cutoff) before mapping to human genes or selecting top 100 available in LINCS
AKI Genes
AKI-associated DEGs were pulled in to compare to the data sets to see if the cystic data sets look more similar to AKI. 157 upregulated and 88 downregulated
aki_upgenes <- c("Adam8", "Adamts1", "Adm", "Akap12", "Akr1b8", "Aldh1a2", "Aloxe3", "Anxa1", "Anxa3", "Apobec3", "Arf2", "Asns", "Atf3", "Bcl3", "Birc3", "Cd14", "Cd44", "Cd68", "Cd9", "Cd93", "Cebpd", "Chrnb1", "Clcf1", "Cldn4", "Cldn7", "Ctgf", "Ctsc", "Cyr61", "Ddx21", "Dusp10", "Dusp5", "Edn1", "Efhd2", "Elf4", "Emp1", "Entpd1", "Epha2", "F2r", "F3", "Fam57a", "Flnc", "Fndc4", "Fosl1", "Fosl2", "Fut2", "Fxyd5", "Gas7", "Gdf15", "Gprc5a", "Havcr1", "Hbegf", "Hilpda", "Hmox1", "Hpcal4", "Igf2bp2", "Il1f6", "Il34", "Irak3", "Itga5", "Klf4", "Klf6", "Klf7", "Krtp28", "Krt20", "Lamc2", "Lif", "Lrp8", "Lrrc32", "Lrrc8c", "Maff", "Map3k6", "Map4k4", "Mapk6", "Mmp3", "Mthfd1l", "Myc", "Myh9", "Nek6", "Ngf", "Nras", "Nt5c1a", "Nucb2", "Oasl1", "Pdgfb", "Pdk4", "Pdlim7", "Plat", "Plaur", "Plin2", "Plk3", "Plp2", "Ppm1j", "Ppp1r14b", "Pprc1", "Procr", "Prrg4", "Ptger4", "Ptgs2", "Ptpn12", "Pvr", "Pxdc1", "Qsox1", "Rab31", "Rad18", "Rap2b", "Rbm3", "Rell1", "Rhbdf2", "Rin1", "Rnd1", "Rnd3", "Rras2", "Rtn4", "Runx1", "Samd4", "Sbno2", "Sema7a", "Serpinb1a", "Serpine1", "Sfn", "Slc16a1", "Slc25a24", "Slc38a2", "Slc7a5", "Slc7a6", "Smad1", "Smox", "Socs3", "Sox9", "Sphk1", "Sprr1a", "Sprr2f", "Sprr2g", "Spry2", "Spsb1", "St3gal1", "Stil", "Syt12", "Taf1d", "Taf4b", "Tcerg1", "Tgfbr1", "Thbd", "Timp1", "Tinagl1", "Tmcc3", "Tmed5", "Tmem173", "Tnc", "Tnfrsf12a", "Tnfrsf1a", "Tpm3", "Trmt61a", "Tubb6", "Txnrd1", "Uchl1", "Vopp1")
aki_downgenes <- c("1700040L02Rik", "Abhd14b", "Acad12", "Acmsd", "Acot11", "Acox2", "Acss2", "Ak4", "Akr1c14", "Akr1c18", "Aldh4a1", "Apeh", "Atp6v1b1", "BC067074", "Bcat1", "Bdh1", "Bhmt2", "Bphl", "Bsnd", "Casr", "Cbs", "Cdkl1", "Ceacam2", "Cmah", "Crot", "Cth", "Cubn", "Cyp2j11", "D3Ertd751e", "Dhdh", "Dnajc28", "Fam107a", "Fggy", "Fmo5", "Fras1", "G6pc", "Galm", "Gatb", "Gatm", "Glyctk", "Gm10804", "Hgd", "Hnmt", "Hsd3b2", "Hykk", "Ift122", "Inpp5j", "Kcnj1", "Keg1", "Klhl3", "Kmo", "Lrrc31", "Lyplal1", "Map2k6", "Mccc1", "Mccc2", "Mep1b", "Mgam", "Nccrp1", "Oxgr1", "Pbld1", "Pde4c", "Pfkm", "Pnkd", "Pter", "Ranbp3l", "Rhcg", "Serpina1d", "Serpina1f", "Sfrp1", "Shmt2", "Slc15a2", "Slc16a9", "Slc22a13", "Slc22ap28", "Slc22a26", "Slc23a1", "Slc25a21", "Slc8a1", "Slco4c1", "Smarca2", "Snx29", "Suox", "Tln2", "Tmem207", "Tmem25", "Tpmt", "Ugt8a")combine all data sets
all_p70 <- read.csv(here("res/deseq2_outputs", "deseq2_p70_fullannotatedres.csv"))
all_p21 <- read.csv(here("res/deseq2_outputs", "deseq2_p21_fullannotatedres.csv"))
all_p28 <- read.csv(here("res/deseq2_outputs", "deseq2_p28_fullannotatedres.csv"))
#at least some lfc and padj > 0.1
all_p70 <- dplyr::filter(all_p70, log2FoldChange > 0.15 | log2FoldChange < -0.15 & padj < 0.1)
all_p21 <- dplyr::filter(all_p21, log2FoldChange > 0.15 | log2FoldChange < -0.15 & padj < 0.1)
all_p28 <- dplyr::filter(all_p28, log2FoldChange > 0.15 | log2FoldChange < -0.15 & padj < 0.1)
combined_fc <- merge(all_p21, all_p28, by = "X")
combined_fc <- merge(all_p70, combined_fc, by = "X")
lfc_aki_up <- combined_fc %>% dplyr::filter( X %in% aki_upgenes) %>% dplyr::mutate(direction = "aki_up") %>% dplyr::select(Genes = X, PrecysticP70 = log2FoldChange, CysticP21 = log2FoldChange.x, CysticP28 = log2FoldChange.y, Direction = direction)#152
lfc_aki_down <- combined_fc %>% dplyr::filter( X %in% aki_downgenes) %>% dplyr::mutate(direction = "aki_down") %>% dplyr::select(Genes = X, PrecysticP70 = log2FoldChange, CysticP21 = log2FoldChange.x, CysticP28 = log2FoldChange.y, Direction = direction)
aki_res <- rbind(lfc_aki_up, lfc_aki_down)AKI Heatmap
##prepare heatmap annotations
aki_direction <- dplyr::select(aki_res, Direction)
aki_direction$Direction <- factor(aki_direction$Direction, levels = c("aki_up", "aki_down"))
aki_res2 <- tibble::column_to_rownames(aki_res, var = "Genes") %>% dplyr::select(PrecysticP70, CysticP21, CysticP28)
samples <- str_sub(colnames(aki_res2), start = 1, end = -4) %>% as.data.frame(row.names = colnames(aki_res2))
colnames(samples) <- "Cystic_Status"
samples$Cystic_Status <- factor(samples$Cystic_Status, levels = c("Precystic", "Cystic"))
#show_col(viridis(8, option = "D")) #selecting colors
my_color = list(Direction = c(aki_down = "#575C6DFF", aki_up = "#FFEA46FF"), Cystic_Status = c(Precystic = "#00204DFF", Cystic = "#A69D75FF"))plot
png(here("res/deseq2_outputs", "heatmap_aki_genes.png"))
#aki_heatmap <-
pheatmap(aki_res2, annotation_colors = my_color, annotation_row = aki_direction, annotation_col = samples, cluster_rows = TRUE, cluster_cols = TRUE, show_rownames = TRUE, color = viridis(n= 10000, alpha = 1, begin = 0, end = 1, option = "viridis"), border_color = "NA", main = "Expression of AKI Genes Across Precystic/Cystic Data", angle_col = "0")## Warning: The input is a data frame, convert it to the matrix.dev.off()## quartz_off_screen
## 2#aki_heatmapwrite.csv(aki_res, file = here("res/deseq2_outputs", "aki_gene_lfc.csv"))Load in (mouse DEG) data
’p70
musdegs_p70_UP <- read.csv(here("res/deseq2_outputs", "deseq2_upgenes_res_p70.csv"))
musdegs_p70_DOWN <- read.csv(here("res/deseq2_outputs", "deseq2_downgenes_res_p70.csv") )’p28
musdegs_p28_UP <- read.csv(here("res/deseq2_outputs", "deseq2_upgenes_res_p28.csv") )
musdegs_p28_DOWN <- read.csv(here("res/deseq2_outputs", "deseq2_downgenes_res_p28.csv") )’p21
musdegs_p21_UP <- read.csv(here("res/deseq2_outputs", "deseq2_upgenes_res_p21.csv") )
musdegs_p21_DOWN <- read.csv(here("res/deseq2_outputs", "deseq2_downgenes_res_p21.csv") )Venn diagram of all DEGs
degsetp70 <- c(musdegs_p70_UP$Mouse.gene.stable.ID, musdegs_p70_DOWN$Mouse.gene.stable.ID)
degsetp21 <- c(musdegs_p21_UP$Mouse.gene.stable.ID, musdegs_p21_DOWN$Mouse.gene.stable.ID)
degsetp28 <- c(musdegs_p28_UP$Mouse.gene.stable.ID, musdegs_p28_DOWN$Mouse.gene.stable.ID)
venn3(degsetp70, degsetp21, degsetp28, categories = c("Precystic P70" , "Cystic P21" , "Cystic P28"), filename = here("res/deseq2_outputs", "degs_venndiagram.png"))## [1] 1FEA p70 precystic ’p70 Mouse
fea_mouse_p70 <- mouseenrich(musdegs_p70_UP$Mouse.gene.stable.ID, musdegs_p70_DOWN$Mouse.gene.stable.ID, 30)
bp_mouse_p70_overall <- barplot(fea_mouse_p70$topcompared) + ggtitle("Pre-Cystic Top Upregulated ")
bp_mouse_p70_overall
mouse_p70_bubble <- bubbleplot(fea_mouse_p70$topcompared)
mouse_p70_bubble
bubbleplot
bubbleplot_p70 <- bubbleplot(dplyr::filter(fea_mouse_p70$topcompared, intersection_size > 6)) + labs(title = "Top Enriched Terms for Pre-Cystic P70 Dataset") + theme(text = element_text(size = 14), plot.title = element_text(size = 24))
bubbleplot(fea_mouse_p70$topcompared) + labs(title = "Top Enriched Terms for Pre-Cystic P70 Dataset") + theme(text = element_text(size = 14), plot.title = element_text(size = 24))
ggsave(bubbleplot_p70, filename = here("res", "fea", "p70_fea_bubbleplot.pdf"), height = 7, width = 14)FEA p21 cystic ’p21 Mouse
fea_mouse_p21 <- mouseenrich(musdegs_p21_UP$Mouse.gene.stable.ID, musdegs_p21_DOWN$Mouse.gene.stable.ID, 30)bubbleplot
fea_mouse_p21$topcompared$term_name <- stringr::str_wrap(fea_mouse_p21$topcompared$term_name, width = 70)
bubbleplot_p21 <- bubbleplot(fea_mouse_p21$topcompared) + labs(title = "Top Enriched Terms for Cystic P21 Dataset") + theme(text = element_text(size = 14), plot.title = element_text(size = 24))
bubbleplot_p21
ggsave(bubbleplot_p21, filename = here("res", "fea", "p21_fea_bubbleplot.pdf"), height = 7, width = 14)FEA p28 cystic ’p28 Mouse
fea_mouse_p28 <- mouseenrich(musdegs_p28_UP$Mouse.gene.stable.ID, musdegs_p28_DOWN$Mouse.gene.stable.ID, 30)bubbleplot
fea_mouse_p28$topcompared$term_name <- stringr::str_wrap(fea_mouse_p28$topcompared$term_name, width = 70)
bubbleplot_p28 <- bubbleplot(fea_mouse_p28$topcompared) + labs(title = "Top Enriched Terms for Cystic P28 Dataset") + theme(text = element_text(size = 14), plot.title = element_text(size = 24))
bubbleplot_p28
ggsave(bubbleplot_p28, filename = here("res", "fea", "p28_fea_bubbleplot.pdf"), height = 7, width = 14)Save gprofiler results
#fea_mouse_p70
fea_mouse_p70_upterms <- apply(fea_mouse_p70$upregulated, 2, as.character)
fea_mouse_p70_downterms <- apply(fea_mouse_p70$downregulated, 2, as.character)
write.csv(rbind(fea_mouse_p70_upterms, fea_mouse_p70_downterms), file = here("res", "fea", "fea_p70degs.csv"))
#feamouse_p21
fea_mouse_p21_upterms <- apply(fea_mouse_p21$upregulated, 2, as.character)
fea_mouse_p21_downterms <- apply(fea_mouse_p21$downregulated, 2, as.character)
write.csv(rbind(fea_mouse_p21_upterms, fea_mouse_p21_downterms), file = here("res", "fea", "fea_p21degs.csv"))
#fea_mouse_p28
fea_mouse_p28_upterms <- apply(fea_mouse_p28$upregulated, 2, as.character)
fea_mouse_p28_downterms <- apply(fea_mouse_p28$downregulated, 2, as.character)
write.csv(rbind(fea_mouse_p28_upterms, fea_mouse_p28_downterms), file = here("res", "fea", "fea_p28degs.csv"))Gene overlaps
Overlap comparing DEGs before human ortholog mapping, using Ensembl ID’s
degs_p28_p21_overlap_up <- merge(musdegs_p28_UP, musdegs_p21_UP, by = "Mouse.gene.stable.ID")
degs_p28_p21_overlap_down <- merge(musdegs_p28_DOWN, musdegs_p21_DOWN, by = "Mouse.gene.stable.ID")
degs_p70_p21_overlap_up <- merge(musdegs_p70_UP, musdegs_p21_UP, by = "Mouse.gene.stable.ID")
degs_p70_p21_overlap_down <- merge(musdegs_p70_DOWN, musdegs_p21_DOWN, by = "Mouse.gene.stable.ID")
degs_p70_p28_overlap_up <- merge(musdegs_p70_UP, musdegs_p28_UP, by = "Mouse.gene.stable.ID")
degs_p70_p28_overlap_down <- merge(musdegs_p70_DOWN, musdegs_p28_DOWN, by = "Mouse.gene.stable.ID")
degs_all_overlap_up <- merge(degs_p28_p21_overlap_up, musdegs_p70_UP, by = "Mouse.gene.stable.ID")
degs_all_overlap_down <- merge(degs_p28_p21_overlap_down, musdegs_p70_DOWN, by = "Mouse.gene.stable.ID")
alldegs_overlap <- rbind(degs_all_overlap_up, degs_all_overlap_down)
alldegs_lfc <- cbind(CysticP28 = alldegs_overlap$log2FoldChange.x, CysticP21 = alldegs_overlap$log2FoldChange.y, PrecysticP70 = alldegs_overlap$log2FoldChange)
row.names(alldegs_lfc) <- alldegs_overlap$X.xPrecyst Unique and Cystic Overlap
Cystic overlap enrichment
cystic_geneoverlap_enrichment <- mouseenrich(degs_p28_p21_overlap_up$Mouse.gene.stable.ID, degs_p28_p21_overlap_down$Mouse.gene.stable.ID, 30)
cystic_geneoverlap_enrichment$topcompared$term_name <- stringr::str_wrap(cystic_geneoverlap_enrichment$topcompared$term_name, width = 70)
cystic_overlap_bubbleplot <- bubbleplot(cystic_geneoverlap_enrichment$topcompared)+ labs(title = "Top Enriched Terms for Overlapped Cystic DEGs") + theme(text = element_text(size = 10), plot.title = element_text(size = 11))
cystic_overlap_bubbleplot
Precystic unique enrichment
#up genes
precystic_uniquegene_enrichment_up <- dplyr::anti_join(musdegs_p70_UP, musdegs_p28_UP, by = "Mouse.gene.stable.ID")
precystic_uniquegene_enrichment_up <- dplyr::anti_join(precystic_uniquegene_enrichment_up, musdegs_p21_UP, by = "Mouse.gene.stable.ID")
#down genes
precystic_uniquegene_enrichment_down <- dplyr::anti_join(musdegs_p70_DOWN, musdegs_p28_DOWN, by = "Mouse.gene.stable.ID")
precystic_uniquegene_enrichment_down <- dplyr::anti_join(precystic_uniquegene_enrichment_down, musdegs_p21_DOWN, by = "Mouse.gene.stable.ID")
precystic_unique_enrichment <- mouseenrich(precystic_uniquegene_enrichment_up$Mouse.gene.stable.ID, precystic_uniquegene_enrichment_down$Mouse.gene.stable.ID, 30)
precystic_unique_enrichment$topcompared$term_name <- stringr::str_wrap(precystic_unique_enrichment$topcompared$term_name, width = 60)
precyst_unique_bubbleplot <- bubbleplot(precystic_unique_enrichment$topcompared)+ labs(title = "Top Enriched Terms for Unique Precystic DEGs") + theme(text = element_text(size = 10), plot.title = element_text(size = 11))
precyst_unique_bubbleplot
save
ggsave(cystic_overlap_bubbleplot, filename = here("res", "fea", "cystic_geneoverlap_fea.pdf"), height = 6, width = 8)
ggsave(precyst_unique_bubbleplot, filename = here("res", "fea", "precystic_geneunique_fea.pdf"), height = 6.5, width =7)Cystic Overlap (’p28 and ’p21)
merge for overlap
fea_cystic_overlap_up <- merge(fea_mouse_p21$upregulated, fea_mouse_p28$upregulated, by = "term_name")
fea_cystic_overlap_down <- merge(fea_mouse_p21$downregulated, fea_mouse_p28$downregulated, by = "term_name")merge for overlap across all datasets
fea_all_overlap_up <- merge(fea_cystic_overlap_up, fea_mouse_p70_upterms, by = "term_name")
#no res for p70 down
#fea_all_overlap_down <- merge(fea_cystic_overlap_down, fea_mouse_p70_downterms, by = "term_name")FEA unique to precystic
#termsinp28 <- fea_mouse_p70$upregulated$term_size[fea_mouse_p70$upregulated]
fea_p70_unique <- dplyr::anti_join(fea_mouse_p70$upregulated, fea_mouse_p21$upregulated, by = "term_name")
fea_p70_unique <- dplyr::anti_join(fea_p70_unique, fea_mouse_p28$upregulated, by = "term_name")
#goterms <- fea$term_id[ fea$source == "GO:BP" & fea$direction == direction]
ggplot(slice_min(fea_p70_unique, order_by = recall, n = 20), aes(x=intersection_size, y=term_name, fill = p_value)) +
geom_bar( stat = "identity") +
scale_fill_continuous(type = "viridis") +
labs(x = "Number of Genes Matched to Term", y = "Functional Terms", title = "Terms Unique to Precystic DEG Enrichment") +
theme_minimal(base_size = 8)
Pathway Similarity
rrvgo package, uses GO terms for input ###### P70
sim_deg_p70 <- go_bp_sim(fea_mouse_p70, "org.Mm.eg.db", threshold = 0.80)## ## preparing gene to GO mapping data...## preparing IC data...#options(ggrepel.max.overlaps = 50)
sp_deg_p70 <- scatterPlot2(sim_deg_p70$simMatrix, sim_deg_p70$reducedTerms, labelSize = 2.5) + labs(title = "Pre-cystic P70 Signature Enrichment Overview")
sp_deg_p70
tm_deg_p70 <- treemapPlot2(sim_deg_p70$reducedTerms, title = "Enrichment for Pre-cystic P70 DEGs")
hm_deg_p70 <- heatmapPlot(sim_deg_p70$simMatrix, sim_deg_p70$reducedTerms, annotateParent=TRUE, annotationLabel="parentTerm", fontsize=6)
P21
sim_deg_p21 <- go_bp_sim(fea_mouse_p21, "org.Mm.eg.db", threshold = 0.80)## preparing gene to GO mapping data...## preparing IC data...sp_deg_p21 <- scatterPlot(sim_deg_p21$simMatrix, sim_deg_p21$reducedTerms, labelSize = 2.5) + labs(title = "Cystic P21 Signature Enrichment Overview")
sp_deg_p21
tm_deg_p21 <- treemapPlot2(sim_deg_p21$reducedTerms, title = "Enrichment for Cystic P21 DEGs")
hm_deg_p21 <- heatmapPlot(sim_deg_p21$simMatrix, sim_deg_p21$reducedTerms, annotateParent=TRUE, annotationLabel="parentTerm", fontsize=6)
P28
sim_deg_p28 <- go_bp_sim(fea_mouse_p28, "org.Mm.eg.db", threshold = 0.80)## preparing gene to GO mapping data...## preparing IC data...sp_deg_p28 <- scatterPlot2(sim_deg_p28$simMatrix, sim_deg_p28$reducedTerms, labelSize = 2.5) + labs(title = "Cystic P21 DEGs Enrichment Overview")
sp_deg_p28
tm_deg_p28 <- treemapPlot2(sim_deg_p28$reducedTerms, title = "Enrichment for Cystic P28 DEGs")
hm_deg_p28 <- heatmapPlot(sim_deg_p28$simMatrix, sim_deg_p28$reducedTerms, annotateParent=TRUE, annotationLabel="parentTerm", fontsize=6)
#ggsave(hm, filename = here("res", "test.pdf"))Term overlaps
Cystic overlap
fea_overlap_p28_p21_up <- merge(fea_mouse_p28$upregulated, fea_mouse_p21$upregulated, by = "term_name")
fea_overlap_p28_p21_down <- merge(fea_mouse_p28$downregulated, fea_mouse_p21$downregulated, by = "term_name")
fea_overlap_p28_p21_up$term_name## [1] "CXCR chemokine receptor binding"
## [2] "fever generation"
## [3] "positive regulation of acute inflammatory response"
## [4] "positive regulation of hair follicle development"
## [5] "positive regulation of monocyte chemotaxis"
## [6] "regulation of acute inflammatory response to antigenic stimulus"
## [7] "regulation of monocyte chemotaxis"fea_overlap_p28_p21_down$term_name## character(0)Precystic overlap
fea_overap_all_up <- merge(fea_mouse_p70$upregulated, fea_overlap_p28_p21_up, by = "term_name")
fea_overap_all_down <- merge(fea_mouse_p70$downregulated, fea_overlap_p28_p21_down, by = "term_name")
fea_overap_all_up$term_name## character(0)fea_overap_all_down$term_name## character(0)Save
fea_overlap_p28_p21_upterms <- apply(fea_overlap_p28_p21_up, 2, as.character)
write.csv(fea_overlap_p28_p21_upterms, file = here("res", "fea", "fea_overlap_p28_p21_upterms.csv"))
fea_overlap_p28_p21_downterms <- apply(fea_overlap_p28_p21_down, 2, as.character)
write.csv(fea_overlap_p28_p21_downterms, file = here("res", "fea", "fea_overlap_p28_p21_downterms.csv"))
fea_overap_all_upterms <- apply(fea_overap_all_up, 2, as.character)
write.csv(fea_overap_all_upterms, file = here("res", "fea", "fea_overap_all_upterms.csv"))
fea_overap_all_downterms <- apply(fea_overap_all_down, 2, as.character)
write.csv(fea_overap_all_downterms, file = here("res", "fea", "fea_overap_all_downterms.csv"))Versions
R.Version()## $platform
## [1] "x86_64-apple-darwin17.0"
##
## $arch
## [1] "x86_64"
##
## $os
## [1] "darwin17.0"
##
## $system
## [1] "x86_64, darwin17.0"
##
## $status
## [1] ""
##
## $major
## [1] "4"
##
## $minor
## [1] "2.0"
##
## $year
## [1] "2022"
##
## $month
## [1] "04"
##
## $day
## [1] "22"
##
## $`svn rev`
## [1] "82229"
##
## $language
## [1] "R"
##
## $version.string
## [1] "R version 4.2.0 (2022-04-22)"
##
## $nickname
## [1] "Vigorous Calisthenics"installed.packages()[names(sessionInfo()$otherPkgs), "Version"]## treemap signatureSearchData ExperimentHub
## "2.4-3" "1.10.0" "2.4.0"
## AnnotationHub BiocFileCache dbplyr
## "3.4.0" "2.4.0" "2.2.1"
## ComplexHeatmap RColorBrewer VennDiagram
## "2.12.1" "1.1-3" "1.7.3"
## futile.logger GOSemSim rrvgo
## "1.4.3" "2.22.0" "1.8.0"
## here dplyr ggplot2
## "1.0.1" "1.0.10" "3.4.0"
## signatureSearch SummarizedExperiment Biobase
## "1.10.0" "1.26.1" "2.56.0"
## GenomicRanges GenomeInfoDb IRanges
## "1.48.0" "1.32.4" "2.30.1"
## S4Vectors BiocGenerics MatrixGenerics
## "0.34.0" "0.42.0" "1.8.1"
## matrixStats Rcpp stringr
## "0.62.0" "1.0.9" "1.4.1"
## viridis viridisLite gprofiler2
## "0.6.2" "0.4.1" "0.2.1"